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Abstract 

On I We get fractional symmetric Fokker - Planck and Einstein - Smoluchowski 



kinetic equations, which describe evolution of the systems influenced by stochas- 
tic forces distributed with stable probability laws. These equations generalize 
On ! known kinetic equations of the Brownian motion theory and contain symmetric 

fractional derivatives over velocity and space, respectively. With the help of these 
j5 ■ equations we study analytically the processes of linear relaxation in a force - free 

case and for linear oscillator. For a weakly damped oscillator we also get kinetic 

. equation for the distribution in slow variables. Linear relaxation processes are 

^ ■ also studied numerically by solving corresponding Langevin equations with the 

O ■ source which is a discrete - time approximation to a white Levy noise. Numerical 

and analytical results agree quantitatively. 
PACS: 05.10 Gg, 05.40. Fb 



1 Introduction 

Studies of evolution of various systems under the influence of external stochastic forces 
constitute content of a large section of statistical physics. It has a great variety of 
applications in physics, chemistry, biology, economy and sociology, see, e.g., |I|, 0, 
0. The most known problem is the description of Brownian motion. In classical 
formulation the Brownian motion manifests itself as unceasing chaotic motion of small 
macroscopic particles in a liquid or a gas. The Brownian motion is explained by 
the presence of unbalanced "pushes" of surrounding atoms and, hence, reveals atomic 
structure of continuous medium, in which the motion occurs. 

The achievements of probability theory serve as a mathematical basis of the kinetic 
theory of the Brownian motion. They are as follows: 

(i) the Central Limit Theorem, which justifies Gaussianity of stochastic force acting 
on Brownian particle; 

(ii) the theory of Markovian stochastic processes; an important assumption used 
at the kinetic description of the Brownian motion is that the behavior of particle at a 
given instant depends only on instantaneous values of physical parameters, but not on 
their previous history; 

(iii) studies of stochastic Gaussian processes, first of all, the work of Bachelier 
(1900), who was the first to study a continuous stochastic Gaussian process with inde- 
pendent increments, and the work of Wiener (1927), who has given a rigorous mathe- 
matical formulation of this process and has studied properties of its sample paths. 

The basic equations of the kinetic theory of Brownian motion are the Fokker - 
Planck equation for the probability density function (PDF) f{x, v, t) in the phase space 
of coordinates and velocities, and the Einstein - Smoluchowski equation for the PDF 
/(x, t) in the real space. The relaxation in the phase space can occur in two steps: the 
first one, "fast" stage, at which the relaxation over velocities occurs, and Maxwellian 
PDF is established, and the second one, "slow" diffusion stage, at which relaxation in 
the real (coordinate) space occurs. If in the system considered the physical conditions 
are such that the two stages of relaxation can be separated, then it is possible to 
pass from the Fokker - Planck equation to the Einstein - Smoluchowski equation and 
describe the system on the diffusion stage with a simpler equation. In more detail the 
transition to the diffusion stage is discussed in the article by Chandrasekhar [^], which 
contains a brilliant presentation of classical theory of the Brownian motion, as well as 
in the monograph containing a modern presentation of the theory of the Brownian 
motion with inclusion of the motion in nonlinear open systems. 

Beginning from the second half of 80-th the term "Levy motion" is becoming widely 
used in statistical physics, in particular, in problems of anomalous diffusion, when char- 
acteristic displacement of diffusive particle grows as t'^ , /x 7^ 1/2 (the case /i = 1/2 
corresponds to the classical Brownian motion). Levy anomalous diffusion appears in 
different areas of physics, namely, turbulence |^ , solid and amorphous states physics 
0, chaotic dynamics [0, plasma physics |§] etc. It is worthwhile to mention also non 
- physical areas, e. g., biology and physiology [H, and the theory of finance |TD[. The 



Levy motion can be considered as a generalization of the Brownian motion. Indeed, 
the mathematical foundation of the generalization are remarkable properties of stable 



probability laws, whose theory has begun from the works of Levy and Khintchine JIT 
From the limit theorems point of view the stable probability laws are the generalization 
of widely used Gaussian law. Namely, stable laws are the limit ones for the probability 
laws of (properly normalized) sums of independent identically distributed (i.i.d.) ran- 



doin variables |T^. Therefore, these laws (as the Gaussian one) occur, when evolution 
of a physical system or the result of an experiment are determined by the sum of a 
large number of identical independent random factors. An important distinction of 
stable PDFs is the power law tail decreasing as |x| "^ , x -^ oo, a is the Levy index, 
< a < 2. Hence, the PDF moments of the order q > a diverge. 

The above mentioned properties of stable PDFs allow one to get a simple intuitive 
basis for anomalous diffusion in the framework of the model of independent random 
"jumps" [^. However, in order to construct a consistent theory of the Levy motion 
it is necessary to get the kinetic equations, which generalize those of the Brownian 
motion, namely, the Fokker-Planck and the Einstein- Smoluchowski equations. It is 
clear also from the very simple arguments that these equations must contain space 
and/or time fractional derivatives. During last two decades a few monographs solely 
devoted to the theory of fractional calculus were appeared, see, e. g., |13[, and an 
extensive treatment of fractional order differential equations applied to heat and mass 
transfer were given ^i\. Different forms of diffusion - like equations with fractional 



derivatives were proposed [T^, |T^, [0, [|18|, ||T9[, pO[. These equations were applied 



in particular, for the description of anomalous diffusion on random fractals 

and in chaotic Hamiltonian systems, for which the orders of space and time fractional 



derivatives are determined by delicate properties of the phase space [0], ||23|. We also 



refer to the paper [Q, where the general description of the fractional relaxation - 
oscillation and fractional diffusion - wave phenomena was provided with adoption a 
mathematical approach to the fractional calculus that is as simple as possible. 

Our paper deals with fractional generalizations of the classical Fokker-Planck and 
Einstein-Smoluchowski equations describing relaxation in the phase and real spaces, 
respectively. We follow classical approach |^, which is used to derive kinetic equations 
for the Brownian motion as well as an approach used in Ref. [ES], where, as far as we 



know, the fractional kinetic equation for the phase space PDF was proposed for the 
first time. Through the paper we restrict ourselves by a one - dimensional coordinate 
space and two - dimensional (coordinate plus velocity) phase space. The generalization 
on two - and three - dimensional real space cases are discussed at the end of the article. 
Besides, we restrict ourselves by symmetric fractional kinetic equations, that is, those, 
which include symmetric fractional derivatives (see below). 

Thus, in our paper we, at first, derive kinetic equations with symmetric fractional 
derivatives, which generalize Fokker - Planck and Einstein - Smoluchowski equations to 
the case of the Levy motion. These equations are called fractional symmetric Fokker - 
Planck equation (FSFPE) and fractional symmetric Einstein - Smoluchowski equation 
(FSESE) , respectively. Secondly, we use these equations for studying classical problems 
of linear relaxation, namely, relaxation in a force - free case and relaxation of a linear 
Levy oscillator. 

It is worthwhile to note that the force - free relaxation in a spatially homogeneous 
case as well as relaxation of a linear oscillator were first studied in Ref. p6|. In this 
paper the equations for characteristic functions were obtained by solving corresponding 
Langevin equations with subsequent averaging the Liouville equation for the phase 
space density. Besides, with the help of fractional kinetic equation for the diffusion 



stage of relaxation the linear oscillator was considered in Ref. ^Tj. With the help 



of FSFPE and FSESE we study in detail both "fast" and "slow" stages of a linear 
relaxation and demonstrate a transition from the first level of description (the use of 
FSFPE) to the second one (the use of FSESE). 

Further, we point attention on a description of two limit cases for the oscillator. 



namely, overdamped oscillator and a weakly damped oscillator. Both cases are very 
important when studying the Levy motion in nonlinear open systems. We propose a 
new kinetic equation for the weakly damped linear oscillator and study its solutions. We 
also solve numerically the Langevin equations, which correspond to fractional kinetic 
equations. We demonstrate numerical results for the linear relaxation problems which 
are solved analytically, and show a good agreement between the results of kinetic theory 
and the results of numerical modelling. 

The paper is organized as follows. In Sec. 2 we derive fractional generalizations of 
the Fokker - Planck and the Einstein - Smoluchowski equations following the approaches 
of 1^, |25|. In Sec. 3 we investigate relaxation in a real space for the force - free Levy 



motion and for the Levy linear oscillator. In Sec. 4 we investigate both relaxation 
problems in the phase space. In Sec. 5 we consider relaxation for the limit cases of 
the Levy oscillator, namely, overdamped and weakly damped oscillators. Finally, the 
conclusions and discussion are presented in Sec. 6. 

2 Fractional Fokker-Planck and fractional Einstein- 
Smoluchowski equations 

It is known that usually the derivation of the Fokker - Planck equation is based on 
the assumption of the finiteness of the second moment of the PDF. Since this assump- 
tion breaks down for the stable PDFs, we find that the classical approach used by 
Chandrasekhar [H] can be adopted for the derivation without using the finiteness of the 



second moment. A similar treatment was undertaken in Ref. 28 , where the discussion 



of Ref. P9 was adopted for the purpose of getting kinetic equation in coordinate space. 



2.1 Fractional Fokker - Planck equation 

Like in the Brownian motion theory, the initial equations in our approach are as follows: 
1) integral equation for the PDF f{x,v,t) of the Markovian stochastic process in 
the phase space: 



f{x,v,t + At) = d{Ax)d{Av)f{x-Ax,v-Av,t)x 

x'i/ix- Ax,v - Av;Ax,Av,At), (2.1) 

where \E'(a:, t;; Aa;, Af , At) is the transition probability, that is, the probability for 
coordinate x to get increment As, and for the velocity v the increment Av during the 
interval At, as well; 



2) the Langevin equations 



dx 

H 
dv 



-UV + F + A{t) , (2.2) 



where u is the friction coefficient, which is assumed to be independent of v, F is the 
regular external force, and A{t) is the fluctuation component of the external force. 

Further in the theory of Brownian motion, by following traditional assumptions [Q , 
one gets expressions for the increments of coordinate and velocity during time interval 



At, which is longer than characteristic time intervals of A{t), though being shorter 
than time intervals during which physical parameters change appreciably: 

Ax = vAt , 

Av = -{uv-F)At + B{At) . (2.3) 

t+At 

Here B{At) = / A(t )dt is a. non-stationary stochastic process, which is assumed to 

t 
be homogeneous Gaussian process with independent increments possessing PDF 

-(i^(At))-^^expf-MM)!) . (2.4) 

The Central Limit Theorem serves as a mathematical justification of this assumption. 
In accordance with stated above, we make a generalization of Chandrasekhar's ap- 
proach with the help of generalization of the Central Limit Theorem to the case of 
i.i.d. variables having infinite variances. Indeed, Levy and Khintchine [^ discovered 
a class of stable probability laws. They are limit ones for the probability laws of nor- 
malized sums of i.i.d. random variables. Each stable law with a characteristic Levy 
index (0 < « < 2) possesses finite moments of orders g, < g < a , but infinite mo- 
ments for higher orders. The Gaussian law is also stable one with a characteristic index 
a= 2, and it possesses moments of all orders. Returning to Eq. (2.4) we note that 
in the theory of stochastic processes the corresponding generalization of homogeneous 
Gaussian process B{t) with independent increments are the Levy stable processes L{t), 
having the characteristic function (we restrict ourselves by symmetric stable laws) ||30| 

WL{k, At)) = (e'^^) = exp(-D \kf At) , (2.5) 

where D > 0, (DAt)^^" is called the scale parameter. At a= 2 one has Gaussian process 
B{t). Above statements justifies expediency of generalization B{At) —>■ L{At) in Eqs. 
(2.3). 

With Eqs. (2.3) and (2.5) the transition probability, see Eq. (2.1) is 

'^{x,v;Ax,Av) =i:{x,v;Av)S{Ax-vAt) , (2.6) 

where S implies the Dirac delta-function, and 

/(] h 
— exp [-ik{Av + uvAt - FAt) -D\k\" At] (2.7) 

— oo 

is the transition probability in the velocity space. 

We insert Eqs. (2.6) and (2.7) into Eq. (2.1), expand the left and the right sides of 
the equation obtained into series of At and set At -^ 0. As the result we get 

^ + v^ = 
dt dx 

d{Av)f{x,v~Av,t) / —exp{—ikAv)x 

— oo 

X [-ikF + ikv{v - Av) + D \k\''] . (2.8) 



Now we turn to the physical space by making an inverse Fourier transform over velocity 
in the right - hand side of Eq. (2.8). We treat each term in the square brackets 
separately. The first and the second terms, being "classical " ones, are transformed 
trivially, leading to 



dv 



and 



"l*''^' 



respectively. We turn to the last term, which can be written as 



/oo/7/<" r^^^ rile .^ 

-oo 2tx j-oo 27r 



(2.9) 



where / (k, /c, t) is the characteristic function, 
f{x,v,t) = 



02 rite rcso rlh ^ 

-00 zvr J-oD /vr 



(2.10) 



We use the symmetric fractional derivative of any order a > 0, which may be 
defined, for a "sufficiently well - behaved" function (p{v),v G R, as the (pseudo - 
differential) operator characterized in its Fourier representation by 



a \v\ 



{v)^-\kf\i){k), v,keK,a>0. 



f2.ir 



Here at the left - hand side we adopt the notation introduced in Ref. [IS 



In order to treat this kind of fractional derivative properly, we remind the definition 



of the left - and right - side Liouville derivatives on the infinite axis |T^ 

1 d r mn 



V^(h(v) 



where < a < 1. For a > 1 



Vl<p{v) 



r(l - a) dv J-00 {x - O'' ' 
r(l - a) dv Jv (^ - xy ' 



f±l^ 



rf' 



Tin — a) dv 



^ POO 



n—a—1 



{vTOd^ 



(2.12) 
(2.13) 

(2.14) 



n = [a] + 1, here the square brackets denote the integer part of the order of the deriva- 
tive. The derivatives (2.12) - (2.14) are characterized in their Fourier representation 
by 



Vl(j){v) ^ (TzA:)"0(fc), a > 0, 



(2.15) 



where 



[T'ik)"' = \k\°' exp ( =F sgnk ) . 



lan 



Thus, the symmetric space - fractional derivative (2.11) can be written as 



d \v\ 



2cos(7ra/2) 



(2.16) 



where a ^ 1,3, ... 

Now we are able to return to Eq. (2.8) and write the kinetic equation for the PDF 
f{x,v,t) in the phase space as 

df df df d d^ 

dt dx dv dv (9|i;|° 

where the last term is defined through Eqs. (2.11) - (2.16). This is a fractional Fokker 

- Planck equation for the Levy motion. At a = 2 this is an ordinary Fokker - Planck 
equation for the Brownian motion. 

Notice that Eq. (2.17) becomes meaningless when a is an odd integer number. 
That is why, for our range of interest < a < 2, the particular case a = 1 must 
be treated separately. However, if one uses the Fourier transform over velocity when 
solving a particular problem, this case is not singled out. 

2.2 Fractional Einstein - Smoluchowski equation 

At the Brownian motion description, along with the relaxation parameter 1/z/, another 
relaxation parameter exists, which characterizes diffusion in a real space. If charac- 
teristic time of this process is much greater than l/u, then it is possible to pass from 
the Fokker - Planck equation for the PDF f{x,v,t) to the Einstein - Smoluchowski 
equation for a simpler PDF f{x,t). 

As in previous sub - section, at the derivation of the fractional symmetric Einstein 

- Smoluchowski equation we follow the reasonings used in the theory of Brownian 
motion. 

Instead of Eq. (2.1), an integral equation in a coordinate space serves as an initial 
one, 

f{x,t + At)=fd{Ax)f{x-Ax,t)'ilj{x-Ax;Ax,At) , (2.18) 

where tplx; Ax, At) is the transition probability, that is, the probability for coordinate 
x to get increment Ax during interval At. 

In the kinetic theory of the Brownian motion the transition to the Einstein - Smolu- 
chowski equation corresponds to neglecting the term dv/dt in the Langevin equation 
(2.2) 0. Thus, instead of two equations (2.2) we have a single Langevin equation, 

^ = - + -m , (2.19) 

dt V V 
and, instead of Eqs. (2.3) we get 

FAt I ,, , , , 

Ax = + -L{At) , (2.20) 

where L{t) is a stable process with symmetric PDF and characteristic function (2.5), 
as before. The transition probability follows from Eq. (2.20), 



FAt\ D 



ilj{x,Ax,At) = / — exp -tkiAx \k\ At 



(2.21) 



We insert Eq. (2.21) into Eq. (2.18), expand left - and right - hand sides of equation 
obtained in Taylor series in t and, then take the limit t ^ 0. As the result we get 



fractional symmetric Einstein-Smoluchowsky equation, 

|.4(£;),£J?^; . ,2.22) 

at ox \v J z/" a |x| 

In the next Sections we give examples of relaxation processes governed by Eqs. (2.22) 
and (2.17). 

3 Solutions to fractional symmetric Einstein-Smo- 
luchowski equation 

In this Section we consider two simple examples of relaxation processes governed by 
FSESE, namely, relaxation in a force - free case and relaxation of Levy linear oscillator. 

3.1 Force - free relaxation. 

Setting F = in Eq. (2.22), we seek for the solution of equation 

df _ D 9" 
dt u°' d \x\ 
with an initial condition 



/, (3.1) 



f{x,t = Q)=5{x-xo). (3.2) 

We pass from Eqs. (3.1) and (3.2) to the equation for the characteristic function /(/t, t), 

(3.3) 



(3.4) 



(3.5) 

Thus, in the force - free case the random process x{t) is a non - stationary Levy stable 
process with independent increments and with the scale parameter 



df 
dt 




with an initial condition 




/(/^,t = 


= 0) = e*'"^'" 


The solution of Eqs. (3.3) and (3.4) is 




/(/t,t) =exp 


< IKXn \k\ t 



{Dt) 



l/a 



In the real space the Levy stable PDFs are expressed in terms of Fox' H functions 



31| . Such representation of all stable PDFs was achieved in Ref. |3^. Mathematical 
details on H functions are presented in Refs. [^, |^. However, in present paper we 
do not touch a real space representation for an arbitrary a. 

Since for the stable PDFs the variance and higher moments of integer order diverge, 
as statistical means characterizing the properties of these processes the moments of 



fractional orders can be used ||2^, |^. In order to guarantee the reahty, they must 



be defined for the modulus of stochastic variable. Therefore, in case of force - free 
relaxation the moments of fractional orders are 



Mx{t;q,a) = {\x — Xo\'^) = / dx\x\'' f {x,t\0) = 

f ((Dt)i/°/i/)5C(g; a), 0<q<a 
\ oo, q > a 



(3.6) 



for < a < 2, whereas 

M.(t;g,2) = ^^C(g;2) (3.7) 

for a = 2 and an arbitrary g, where 

C{q;a)= dx2\x2\'^ -r^exp{-ixiX2-\xi\°') . 

The coefficient C{q; a) can be evaluated with the use of generahzed function theory 



C(q;a) = — sin (—)T(l + q)T(l --) , 0<q<a . (3.8) 

nq V 2 / V a J 

Equations (3.6) - (3.8) have a direct physical consequence for description of an 
anomalous diffusion. Indeed, for the ordinary Brownian motion the characteristic dis- 
placement of a particle may be written through the second moment as 

Ax{t) = M]/2(t; 2, 2) oc t^/^ (3.9) 

One may note from Eqs. (3.7), (3.8), that for the normal diffusion 

M2/'^(t;g,2)oct^/' (3.10) 

at any q and, thus any order of the moment can serve as a measure of a normal diffusion 
rate: 

Ax{t) ^ Ml'\t- g, 2) oc t^'\ (3.11) 

if one is interested in time - dependence of the characteristic displacement, but not in 
the value of the prefactor. We remind that usually just the time - dependence, but not 
the prefactor, serves as an indicator of normal or anomalous diffusion p. In analogy 
with Eq. (3.11) it follows from Eqs. (3.6), (3.8) that the quantity My«(t;g,a) at 
< a < 2 and any q < a can serve as a measure of an anomalous diffusion rate: 

Ax(t) ^ My'?(t; g, a) oc t^/", < g < a< 2. (3.12) 

Here we have the case of a fast anomalous diffusion, or superdiffusion. 

In order to illustrate how C{q\ a) influences absolute value of characteristic dis- 
placement Ax(t), in Fig.l we show the g - th order root of C(g; a) as a function of 
g at different values of a. It is seen that C^/'^ weakly depends on g and is near to 1, 
if one does not choose g close to a. This is especially clearly seen for a > 1. This 
figure illustrates that the value of prefactor, at sufficiently small g, practically does not 
influence the anomalous diffusion rate. 

It is expedient to give two remarks on the above. 

The flrst one is related to applicability of the diffusion - like equation (3.1) at the 
Levy indexes less than unity. Indeed, at < a < 1 the characteristic displacement 



(3.12) grows faster than in the baUistic regime, that is, in case of a free particle motion 
with a finite velocity. This property is often proclaimed as unphysical, and the conclu- 
sion follows about inapplicability of the diffusion - like equation (3.1) at a < 1, see, e. 
g., [^. However, as a counterargument, one could mention the law of relative diffu- 
sion (sometimes called as Richardson law) in the field of locally isotropic turbulence, 
according to which the square of displacement grows as t^, that is, faster than in the 
ballistic regime. 

The second remark is related to the divergence of the second moment M^(t;2,a;) 
for the hyper diffusion. In Ref . p^ , in order to extract the scaling form for the second 



moment, the "walker" is enclosed in an "imaginary growing box". This formal proce- 
dure leads to the diffusion scaling t^/" for the variance, which is consistent with Eq. 
(3.12), and it was implemented numerically. However, it seems that more physically 
relevant procedure must take into account the finite velocity of a diffusive particle. 
This problem is beyond the scope of our paper. We only mention a recent extensive 
discussion on this interesting and important theme [P7|]. 

The results of numerical simulation of the Langevin equation (2.19) are shown in 
Figs. 2 and 3. Here and below the stochastic source A{t) is represented in numerical 
simulation as a discrete approximation of a "white Levy noise" , that is, as a stationary 
consequence of independent identically distributed variables having symmetric stable 
PDF with the Levy index a and the scale parameter equal 1. To obtain the sequence 



we use the generator, which is based on the method of inversion |^ along with the 



Gnedenko limit theorem |12|. This generator was described in our recent papers [39 
40| in more detail. Time interval between the terms of the sequence is equal unity. In 
a force - free problem we estimate numerically the moments Mx{t;q,a) by averaging 
over realizations of x{t). The total number of realization is equal 500, each of length 
512. 

In Fig. 2 we show Mx{t] q, 1) versus t at different g in a log - log scale. At g < 
a = 1 the dependence is well fitted by straight line whose slope allows one to define 
the diffusion exponent. At g > a theoretical value of the moment is infinite, and 
in numerical simulation the moment strongly fluctuates, thus it is unable to get the 
diffusion exponent. 

In Fig. 3 we show the exponent /i in the relation 

M^{t;q,a) oct^" (3.13) 

versus the Levy index a of the discrete approximation of a white noise. The order q 
of the moment is equal 0.25, which is smaller then the smallest value a = 1 used in 
numerical simulation. Theoretical dependence 0.25/a is shown by dotted line. The 
values of /x obtained in simulations are shown by black points. A good agreement 
between the theory and numerical simulation is obvious. 

3.2 Relaxation of linear Levy oscillator 

Setting F = —cu'^x in Eq. (2.22), we seek for the solution of equation 



with an initial condition 



f{x,t = 0)=5{x~Xo). (3.15) 



The equation for the characteristic function is 



ITT = 1^^ T K« /, (3-16) 

Ot V OK Z/" 



and an initial condition is 



/ (k, t = 0) = e*'^^" (3.17) 

The solution of Eqs. (3.16) and (3.17) is obtained by the method of characteristics: 

/(fi;,t) = exp<^mxoe"'^* - L)osc(t) |«r ^ , (3.18) 

where 

D ( ac.2 



^osc it) = -^-—^ 1 - e-— * . (3.19) 

This result was obtained in [^. 

We see from Eqs. (3.18) and (3.19) that the relaxation of oscillator, contrary to 
the force - free case, is not a Levy stable process with independent increments. It can 
be named as "stable - like" or "Levy - like " process, since there exists \k\" in the 
exponent of the characteristic function, however, the scale parameter for the oscillator, 
{Dose (t)) , does not follow as t^/", which is a manifestation of the Levy stable process 
with independent increments, see Eq. (2.5). The process x{t) behaves as a Levy stable 
one only asymptotically at small times, 

t«T^ = ^ , (3.20) 

when the exponent in Eq. (3.19) can be expanded into power series, and we get the 
force - free result with accuracy up to linear in t terms inclusively. On the other hand, 
at t >> Ta; the process x(t) becomes asymptotically stationary process with the stable 
PDF which does not depend on t and with the Levy index a and the scale parameter 

Dli:{t = oo)=[^^Y" . (3.21) 

The PDF f{x,t) is expressed in terms of elementary functions in two cases: 
(i) a = 2 (Brownian oscillator). 



f{x,t) = exp 



2./ \2 ^ 

1 I l^ " ^oe '^ '" 



and 

(ii) a = 1 (Cauchy oscillator). 



AnDosc{t) 1 4:Dosc{t) 

D 



(3.22) 



'2uj^tlv 



f{x,t) = Doscit)/^ 

DUt) + {x- xoe--^^/^f 
Dosc{t) = 4 (l - ^'^'"'^ 



Equations (3.22) and (3.23) are the Gaussian and Cauchy PDFs, respectively, whose 
scale parameters does not depend as t^/^,t, respectively. We note also that only for 
the Brownian oscillator stationary solution has Boltzmann form, 



/,t(x;a = 2) = y'^exp|-^x2| . (3.24) 

Below we return to the problems connected with stationary solutions of fractional 
kinetic equations. 

It is convenient to define fractional moments after subtraction of a stochastic quan- 
tity X its regular part containing initial condition, that is, Xoexp(— u;^t/z/). Thus, the 
moments are 



M^{t;q,a) = { x - Xoe'''''/'' )= dx {xl" f {x , t\0) 

= Dli^it)Ciq;a), (3.25) 

where C {q; a)is the same as in (3.8). 

Numerical simulations of a linear oscillator relaxation include solution of the Langevin 
equation (2.19) with an external force F = —u'^x and calculation of the g-th order 
moments. In Fig. 4 we present the results obtained for various values of u by averag- 
ing over 300 realizations, each of length 4096. The Levy index a is equal 1, and the 
order of the moment is 0.25. The values obtained in numerical simulation are depicted 
by black points whereas the solid line indicates the values estimated from Eq. (3.24) 
by solving FSESE. The vertical mark indicates the value r^, after which the random 
process x{t) becomes stationary one. At t > r^; the moment tends to a constant value 
which is estimated from Eqs. (3.21) and (3.25). Numerical results demonstrate a good 
agreement with theoretical estimates on both non - stationary and stationary stages of 
evolution as well. 

4 Solutions to fractional symmetric Fokker - Planck 
equation 

In this Section, as in previous one, we consider the same examples of relaxation pro- 
cesses, but governed by FSFPE. 

4.1 Force - free relaxation. 

Setting F = in Eq. (2.17), we seek for the solution of equation 

df df d (9" 

with an initial condition 

/ (x, t;, t = 0) = 5 (a; - Xo) 6iv - Vo). (4.2) 

For clarity, it is expedient to consider space - homogeneous relaxation at first. 



4.1.1 Space - homogeneous force - free relaxation. 

We seek for the solution f{v,t) of an equation 

df d 9" 

lir = '^jr^vf)+D—^f (4.3) 

at ov d\v\ 

with an initial condition 

fiv,t = 0)=5iv-vo) (4.4) 

We pass to the characteristic function f{k,t), 

/oo ijU ^ 
-f{k,t)e-^'^ , (4.5) 

-oo ZTC 

which obeys an equation 

|+„,|._fli,r/ (4.6) 

with an initial condition 

/(fc,t = 0)=e''="o (4.7) 

The solution of Eqs. (4.6) and (4.7) is obtained by the method of characteristics: 

f{k,t) = exp [ikvoe-"' - Ikl^Dfjit)} (4.8) 

where 



Z)M(t) = _(l_e— *) , (4.9) 

"ff" means "force - free" case. The space homogeneous relaxation in a force - free 
case was first considered in [^, where Eq. (4.6) was obtained and solved with an 
initial condition (4.7). The process of relaxation is not a Levy stable processes with 
independent increments, but, following terminology used in previous Section, can be 
named stable - like, or Levy - like process, since its characteristic function (4.8) contains 
\k\°' in the exponent, but DpAt) is not a linear function of t, hence scale parameter 

Drf(t)j does not scale as t '". The stable process with independent increments 
arises at small times, 

t«T^ = 1/az/ , (4.10) 

when the exponent in (4.9) can be expanded into power series. With the accuracy up 
to linear in t terms inclusively, we get the Levy stable process. On the other hand, 
at t >> Tj, the stochastic process v{t) becomes asymptotically stationary process with 
the stable PDF independent of t and with the Levy index a and the scale parameter 

K'(« = -)) = (^)"° . (4.11) 

In elementary functions the explicit expressions for f{v,t\vQ) can be obtained in two 
cases: 

(i) a = 2 (Brownian space - homogeneous relaxation), for which (e.g., @]) 



f{v,t) 



^(1-e-^-) 






and (ii) a = 1 (Cauchy space - homogeneous relaxation), for which 

p/ X Du 1 — exp(— i/t) /. ,„x 

vr ^^[v — Vo exp[—ut))'^ + iJ^^l — exp(— z/rj)^ 

The stationary Maxwell PDF is obtained for a = 2 only: 

Here it seems expedient to discuss the problems related to stationary solutions of 
fractional kinetic equations. In the classical theory of Brownian motion equilibrium 
Maxwell PDF over velocity is reached at t — > cxd. It is characterized by the temperature 
T of surrounding medium. There exists relation between parameter D of the PDF of 
the random source in the Langevin equations, see Eqs. (2.2) - (2.4), and the friction 
coefficient z/ : 

D = ^ , 4.15 

m 

where m is a mass of Brownian particle, ks is the Boltzmann constant. The tempera- 
ture T is a measure of a mean kinetic energy of the Brownian particle: 

m(t>^) kBT 
{Ek^n) = ^ = ^- • (4.16) 

Equations (4.15) and (4.16) are examples of fluctuation - dissipation relations. For this 
case the source A{t) in the Langevin equation is called a source of internal fluctuations. 
Relations (4.15) and (4.16) may be not fulfilled, as it takes place, e.g., in auto-oscillation 
systems. In such a case one says that A{t) is a source of external (relatively to the 
system considered) fiuctuations in Eq. (2.2). However, Maxwell - Boltzmann expo- 
nential form of stationary solutions retains 0. As to the Levy motion, the fiuctuation 
- dissipation relations can not be fulfilled, at least, because of the infinity of square 
velocity: (f ^) = oo for < a < 2. Therefore, we can only say about A(t) as about 
the source of external fiuctuations. Moreover, it stems from the example considered in 
this subsection, as well as from the example of linear oscillator considered above, that 
the stationary solutions do not possess Maxwell - Boltzmann form but instead more 
general form of stable distributions. In present there is no theory of equilibrium state 
basing on stable PDFs. Perhaps, the achievement will be done with the help of Tsallis' 



statistics and its generalizations, see recent review ^T[ and references therein. 



We also write the g - th order fractional moment of the velocity. 



M,{t;q,a) = {v-Voe-''' ) = [Dy^{t)) C{q;a) , (4.17) 

where DrHt) is determined by Eq. (4.9). This formula is compared with the results 
of numerical simulations at the end of Subsec.4. 

4.1.2 Space - inhomogeneous force - free relaxation. 

We turn to the force - free relaxation in general case, which is governed by Eqs. (4.1) 
and (4.2). We pass to the characteristic function /(k, A;,t), 

dx 



/<^ fix ^ 

—e-^^^-^'^f{^,k,t) , (4.18) 



which obeys an equation 



with an initial condition 



| + (.,-.)| = _Z,|.|"/. (4.19) 



f{K, k,t = 0) = e^'^^o+^'^^'o . (4.20) 

The solution of Eqs. (4.17), (4.18) is obtained by the method of characteristics 



/(k, /c,t) = exp{iK{xo -\ ) + ivQe {k ) — 



-D / dr 
Jo 



- + {k- - e 



} (4.21) 



Equation (4.19) admits an analytical inverse Fourier transform in case of Brownian 
relaxation, 



(4.22) 
where 



27rAV2--^l 2A^"'' 


— z,/ 


V = v-voe-"', 




A = X — Xq (1- 

V 


-e"'^* 


A = ah-h^. 





This result was presented in 0]. 

For an arbitrary a, < a < 2, the analysis of Eq. (4.19) becomes rather compli- 
cated. However, since we already have an information about velocity relaxation, we 
study evolution of a simpler function. 



f{x,t\xo,vo)= dvf{x,v,t\xo,vo) , (4.23) 

whose characteristic function is 

f{K,k = 0,t\xo,vo) =expliKXo + iK—{l-e-''^)-Dff\t)\K\'^\ , (4.24) 

where 

Dff\t) = ^ J^ril - e-n'' . (4.25) 

As well as in previous Sections, we may call the random coordinate in the process of 
space - inhomogeneous relaxation as a stable - like process. The integral in Eq. (4.25) 
can be evaluated in terms of elementary functions in two cases, 
(i) a = 2 (Brownian motion), for which 



and (ii) a = 1, for which 



Dff\t) = ^{t 



D 1 - e-"^ 



At large times t » t^, where r^, is defined by Eq. (4.10), we get 



iW^ 



Dt 



Df>{t ^oo)^- (4.26) 

and the characteristic function (4.24) coincides with the solution of the Einstein - 
Smoluchowski equation in a force - free case (with xq substituted by Xq + Vq/u), see 
Eq. (3.5). Thus, in the limit of large times the random process x{t) becomes a - 
stable process with independent increments and with the Levy index a and the scale 
parameter [DtY^"' /v. One can see that the space - inhomogeneous relaxation occurs 
in two stages, namely, the fast kinetic stage, at which stationary stable velocity PDF 
is established after the time period r^, and slow diffusion stage, whose characteristic 
relaxation time Tx can be defined, if one introduces an external scale L of the system 
considered: 

r.^^ • (4.27) 

For large enough systems t^ » r^. 

For the coordinate we also write the q - th order fractional moment which is esti- 
mated in numerical simulations below: 

Mx{t- q-a) = (jx-Xo- ^(1 - 6"^^*)^ = {D%\t)y/-C{q; a) . (4.28) 

Numerical simulation of a force - free relaxation process described by FSFPE in- 
cludes the solution of the Langevin equations (2.2) with F = and estimation of the 
moments M^it; q, a) and M^(t; g, a). The results are shown in Figs. 5-7. 

Figures 5 and 6 have an illustrative purpose. In Fig. 5 we show typical velocity 
trajectories (at the left) and coordinate trajectories (at the right) for various values 
of the Levy index. In the figure z/ = 0.03 and, thus the velocity relaxation time r„ 
is equal 20 for a = 1.7, 26 for a = 1.3 and 33 for a = 1.0. On the main part of the 
realizations presented the process v{t) is stationary. In the velocity realizations large 
outliers are clearly seen, which appear due to power asymptotics of the stable PDFs of 
the velocity. With the Levy index decreasing ( from top to bottom) the asymptotics 
become flatter which leads to the growth of outlier amplitudes. Large outliers of the 
velocity, in turn, leads to large jumps on the trajectories x{t), that is. Levy flights, see 
illustrations on the right. 

In Fig. 6 we depict typical trajectories of the velocity (at the left) and of the 
coordinate (at the right) for various values of the friction coefficient z/. The Levy index 
is 1.3. The velocity relaxation time r^ is equal 8, 26, and 77 for the top, middle and 
bottom figures, respectively. For these values the process v{t) is stationary almost on 
the whole length of realization. In the figure the difference between trajectories with 
different r^ are clearly seen: with r^, increasing the relaxation of the velocity outliers 
reduced, whereas trajectories of x(t) become less cutted up and more smoothed. 

In Fig. 7 we depict the velocity moments My (at the left) and coordinate moments 
Mx (at the right) versus time at different values of friction coefficient z/. The moments 
are obtained by averaging over 50 realizations, each of length 4096. The momentum 
exponent q is equal 0.25 for the velocity and the coordinate as well, the Levy index is 
equal 1.3. The moments estimated by numerical solution of the Langevin equations are 
shown by black points, whereas the moments estimated with Eqs. (4.17) and (4.27) 
are shown by solid lines. The vertical mark indicates the relaxation time r^,. At the 
intervals greater than Ty the random process v{t) becomes stationary and the velocity 



moment tends to the constant value D /av, see Eq. (4.9). At the same time it follows 
from right figures that the process x{t) remains non - stationary one, and the moment 
of the coordinate tends to the linear (in a log - log scale) asymptotics, which has a 
slope q/a and shown by dotted line in the right figures. From Fig. 7 we can make a 
conclusion about agreement between the theoretical results obtained for a force - free 
relaxation by solving FSFPE and the results obtained by numerical solution of the 
Langevin equations. 



4.2 Relaxation of linear Levy oscillator 

Setting F = —uP'x in Eq. (2.17), we seek for the solution of equation 

at ox ov ov o\v\ 



(4.29) 



with an initial condition 

f{x,v,t = 0) = S{x- xo) S{v- vq) . (4.30) 

Making Fourier - transform, we get for the characteristic function an equation 



-— + [uk — k)—- + uj k—- = —D \k\ J 



dk 



dn 



(4.31) 



dt 
with an initial condition 

/(fi;,A;,t = 0) = e*'^^«+*^^V (4.32) 

The solution of Eqs. (4.31) and (4.32) is obtained by the method of characteristics: 
/(/t, k,t) = exp{iiu^xoi-k,e-^^' + k2e~^'') + ivo{-p2kie~'P'' + Pikie-"'')- 



D / dr 
Jo 



Pik^e-P'^ - P2k,e-P'^ } 



(4.33) 



where 



u 
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Pi -P2 



UJ" 



We may introduce 



and rewrite Eq. (4.33) in terms of hyperbolic functions as 



/(/t, /c, t) = exp{iciJ^a;oe "^'"^ 
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(4.34) 



This expression is valid for v^ — Aui'^ > . If z/^ — 4cj^ < ,then one introduces 



Ui 



u^- 



v" 



and makes the following changes in Eq. (4.34): 

cosa;it , 



ch— 



1 i^it 
— sli — 
z/i 2 

1 
-shvit - 



1 



2uoi 
1 



sincjit 



In aperiodic case, v = 4a;^ ,evidently, 



2ij- 



sin 2ujit - 



c/i 






— sh 
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■t. 



For the particular case of Brownian oscillator, a = 2, one can get the Gaussian 
PDF in the phase space by using an inverse Fourier transform of Eq. (4.34). Further, 
by using the first and the second derivatives of the characteristic function at k = A; = 0, 
one can get the means and the variances of the velocity and coordinate as well. For the 
Brownian oscillator these formula are presented in [^. The expressions for the means 
are also valid for all a greater than unity. 

Now we turn to the more complicated fractional moments: 






(4.35) 



where 



Ditit) = D f dre---^^/' 
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(4.36) 



and C(g; a) is the same as in previous Sections. Equations (4.35) and (4.36) are 
compared with the results of numerical simulation at the end of Sec. 5. 

Equation (4.34), in principle, allows one to study stationary solution, which is 

defined by 

r°o dn r°° dk 
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(4.37) 



(4.38) 



A simple analytic expression can be obtained for the stationary solution in the case 

a = 2 only: 

r n /..2 

.2 



fst{K, k) = exp <^ -— \^ + k 
2v Xw^ 



and 



fst{x,v) 



UUJ 



V 



2^2^ 



2vrD^"Pr2D^" +^^ 



However, from Eqs. (4.37) and (4.38) one can get some conclusions for simpler sta- 
tionary PDFs, namely, 

dvfst{x, v) = / —e-'^'^fstin, A; = 0) , 



and 



—e-'''^UK = 0,k) 

-oo zvr 



Both stationary PDFs are stable ones with the Levy index a and scale parameters, 
which are expressed as integrals over r, see Eq. (4.36). 

Though linear oscillator, as we see, admits a strict solution, the general formulae 
are not easy to analyze analytically. Therefore, it is instructive to consider two limit 
cases, namely, overdamped oscillator and a weakly damped oscillator. Both cases also 
are very important in problems related to a non - linear oscillator. 

5 Limit cases of Levy oscillator 

5.1 Overdamped Levy oscillator, u/u <C 1. 

Let us consider the relaxation of the moments in case of an overdamped oscillator. 
First we turn to the velocity relaxation. It follows from Eqs. (4.35) and (4.36), that we 
can restrict ourselves by zero - order approximation in u/u. In this case we get from 
Eq.(4.36) 

Dl^ilit) ^ ^ {l - e--'^^) , (5.1) 

which is, of course, coincides with DrAt) for the force - free case, and the velocity 
relaxation time r^, for an overdamped oscillator is given by Eq. (4.10). The conclusion 
is that the velocity relaxation for the overdamped oscillator in the main order in small 
parameter uj/v is the same as in the force - free case. 

Now we consider space relaxation, which differs from the force - free case. We get 
from Eq.(4.36) in the first order in uj/v, 

Ditit) ^ ^ 1^ dr {e-'^/'^ - e-^^y (5.2) 

We note that at cj = we get Eq. (4.25) describing relaxation in a force - free case. It 
follows from Eq. (5.2) that at 

t > r^, = l/au 

the second term in the brackets contributes negligibly, and 

^i'slit » r.) ^ ^^, {l - e—'^/'^) , (5.3) 

which coincides with the result obtained in the framework of FSESE, see Eq. (3.19). 
Thus we conclude, that at time intervals greater than the velocity relaxation r^,, the 
overdamped oscillator can be described with the help of a more simple kinetic equa- 
tion, namely, FSESE. For an overdamped oscillator the relaxation process occurs in 
two stages: fast velocity relaxation stage, at which stationary stable velocity PDF is 
established during time interval r^,, and slow diffusion stage, at which during time 
Tx ,see Eq. (3.20), the stable PDF is established in a real space. 



5.2 Weakly damped Levy oscillator, uj/u ^ 1. 

For this case in the theory of Brownian oscillator there exists the method of simplifying 
kinetic description . It is based on the method of slowly varying amplitudes, or the 
van-der-Pol method, which is frequently used, e.g., in radiophysics ||1^, [|^. In this 
approach the solution of the Langevin equations 

dx 



is looking in the form 



— = -uj\x -iyv + A{t) (5.4) 

(Jjv 

V 

X = X COS c^o^ H sm uot , 

V = vcos ujot — uqx sin uot , (5.5) 

where x and v are slowly varying (during the period 2tt/uj ) amplitudes. The choice of 
the solution (5.5) is equivalent to the condition 

dx sin cut dv ^ ,_ „, 

cosLot— — I — = . (5.6) 

dt uj dt 

We insert (5.5) into (5.4) and, after averaging over the period, get the Langevin equa- 
tions for X and v : 

dx V _ , . . 

di + r = ^ • <5-^' 

where 

1 /■* 
A~{t) = / dilA{i!)m\uji! , 

A~{t) = —j dt'A{t') cos ujt' . (5.8) 

2n Jt-2iT/uj 

It follows from (5.7) that the Langevin sources A~{t) and A~{t) does not contain "fast" 
time 2tt/uj. Therefore it stems from Eq. (5.8) that A{t) can be represented as 

A{t) ^ a(t) cos ujt - h{t) sin ojt , (5.9) 

where a{t) and h{t) are random stationary functions which are related to A~(t) and 
A~{t) as 

a{t) = 2A~{t) , 

h{t) = -2ujA~{t) . (5.10) 

Equations (5.9) and (5.10) have the following meaning ^^. According to Eqs. (5.7) 
and (5.8) the random force influences oscillator by means of slowly varying components 
A~{t) and A~{t) (or a{t) and h{t)) only. Therefore, if one considers the random influ- 
ence on the weakly damped oscillator, the main components of the random force are 
singled out by Eq. (5.9). Further, if A(t) is a stationary Gaussian process, then after 
constructing expression for the correlation function of this process, one can convince 
himself that one - point PDFs of A{t), a{t) and b{t) coincide ^^. We assume that the 
conclusion about identical PDFs of A(t), a(t) and b{t) is also valid for (symmetrical) 



stable PDFs notwithstanding the fact that the proof of this statement is not so trivial 
as in the Gaussian case. It follows from this coincidence that the PDFs of the processes 



dt'Ait' 

i-t+At 

La{At) = / dt'a{t' 



t+At 



Lfe(At) = / dt'b{t') (5.11) 

also coincide, and thus, with taking Eq.(2.5), the PDFs of La , Lf, are 

/oo (^h 
— exp{-tkLa,b-D\k\''At) . (5.12) 
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We also define the processes 

/•t+At 

Lz{At) = j^ dt'A~{t') , 

/•t+At 

L~(At) = / dt'A~{t') , (5.13) 



which, according to (5.10) are related to La-, Lb as 

L„(At) = 2L^(At) , 

LbiAt) = -2iuL~{At) . (5.14) 

Now, with the help of Eqs. (5.12) and (5.14) we are able to get the characteristic 
functions w{L~) and w{L~) as well as their PDFs w{L~) and w{L~) : 

wiL^)=expi-D~\krAt) , 

w{L~)=expi-D~\krAt) , (5.15) 

where 

(2c^)" ' ^ 2" 

Equation for f{x,v,t) is derived in the manner analogous to FSESE. An initial 
equation is 

f{x,v,t + At) = f f d{Ax)d{Av)f{x- Ax,v- Av,t) X 

x'^ix- Ax,v- Av;Ax,Av,At) , (5.17) 

where ^ is the transition probability. For the increments Ax, At; we get from Eqs. 

(5.7) 

Ax + '^xAt = L~{At) , 

A^ + ^^At = L~(At) , (5.18) 

where the PDFs for L~ and L:^ are given by Eqs. (5.15) and (5.16). Now we construct 
\E'. From the structure of Eqs. (5.18) it follows that 

^!{x- Ax,v- Av]Ax,Av,At) = ^!~{x- Ax]Ax,At)'^~{v- Av]Av,At) , (5.19) 



where 



^-(2; Ax, At) = [ — exp ( -m [Ax + -xAt) - D~ |/t|" At 



— exp -m Av + -vAt ) - /^~ |A;|" At 



(5.20) 



Insert Eqs.(5.19) and (5.20) into Eq. (5.17), expand into power series in At and tend 
At to zero. As the resuh we get 



dt 



/OO Wtj- /"OO w/p 
— / — exp (—inAx — ikAv) x 
-OO 2% J -OO 2tt 



iK^{x - Ax) + ik^{v - Ad) + D~ |/t|" + D~\k\' 



(5.2i; 



Transforming the terms in the right hand side as it is described in Sec. 2, we arrive at 
the differential equation for f{x,v,t) : 



We find the solution of Eq. (5.22) with an initial condition 

f{x, v,t = 0) = 6{x — xo)6{v — vq) 
We pass to the characteristic function /(k, k,t), 

f{x,v,t)= -- --exp{-iKX-ikv)f{K,k,t) 

j-oo Zir J-oo 111 

which obeys an equation 

df i^ df u df nil"? n 171" f 

9^ = -2^9;: -2^9^-''^""" /-^d^i / , 

with an initial condition 

/(/t, A;, t = 0) = e^'^^o+^^^o . 

The solution of Eqs. (5.25) and (5.26) is 

/(/t, k, t) = exp{iKXoe~'^*''^ + ikvoe~'^*''^ 
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We also get the fractional moments: 



where 
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(5.22) 
(5.23) 

(5.24) 

(5.25) 
(5.26) 



(5.27) 



(5.28) 



(5.29) 



It follows from the last equations that, opposite to the case of overdamped Levy oscil- 
lator, the relaxation times for the weakly damped oscillator are equal, 

2 

T"i, = T"x- = T" = — , (5.30) 

av 

and, thus for a weakly damped oscillator it is impossible to distinguish kinetic and 
diffusion stages of relaxation. At time intervals greater than r the random processes 
x{t) and v{t) becomes stationary ones with stable PDFs. The characteristic function 
of the stationary state is defined from Eq. (5.27): 

2D~ . ,„ 2Dz 



U'^,k)=exp( ^l^r '-\k 



and the PDF retains Maxwell - Boltzmann form for a = 2 only. 

Numerical simulations of a linear oscillator relaxation include solution the Langevin 
equations (2.2) with an external force F = —uj'^x and subsequent calculation of the 
velocity and coordinate moments. The results are shown in Figs. 8 - 13. 

Figures 8-10 has an illustrative character. In Fig. 8 the typical trajectories 
are depicted of the velocity (at the top) and of the coordinate (at the bottom) for 
the overdamped oscillator (at the left) and for the weakly damped oscillator (at the 
right), respectively. The frequency value is equal 0.003 and 0.3 for the overdamped 
and weakly damped oscillators, respectively. The friction coefficient is equal 0.03 and 
the Levy index is equal 1.3. In the figures the trajectories are shown which have a 
single large outlier. It allows us to demonstrate visually the difference in behaviors 
of two kinds of oscillator: the relaxation process for an overdamped oscillator (at the 
left) resembles relaxation in a force - free case and is radically different from rapidly 
oscillating behaviour of the velocity and coordinate of a weakly damped oscillator (at 
the right). 

In Fig. 9 trajectories on the phase plane {x, v) are depicted for (a) weakly damped 
oscillator, u = 0.06, and (b) overdamped oscillator, u = 0.009, respectively. The 
friction coefficient u is 0.02, and the Levy index a is 1.2. In simulations a single 
trajectory having the length 12.000 is used. In both cases the trajectory starts from 
the point (0, 0), and, then is "thrown away" from this point by Levy "pushes" produced 
by external source in the Langevin equation (2.2). For the weakly damped oscillator 
the whole picture is a set of spirals gathering in a focus at (0, 0). For the overdamped 
oscillator the phase trajectories are a set of curves gathering in a node at (0, 0) without 
circumvolution around it. In Fig. 10 the trajectories are depicted on the plane {x, E) 
where E is the energy of oscillator, E = v'^/2 + {ux)'^/2, for the two Levy indexes 
a = 1.9 (above) and a = 1.6 (below). In simulations we use the following parameters: 
u = 0.003, u = 0.002, the length of trajectory is 4096. In the top figure the large jumps 
are almost absent, and the trajectory resembles that of the Brownian oscillator. In the 
bottom figure a single large jump, that is "Levy fiight" occurs, which is due to slowly 
decreasing power law asymptotics of stable distribution. 

In Fig. 11 we show the velocity moments M^ (in the top) and coordinate moments 
Mx (in the bottom) versus t for the overdamped oscillator (at the left) and for the 
weakly damped oscillator (at the right), correspondingly. The oscillation frequencies u 
are equal 0.01 and 0.1 for the overdamped and the weakly damped oscillators, corre- 
spondingly, the friction coefficient is equal 0.03. The order of the moment is 0.5, and 
the Levy index is 1.3. The moments obtained by numerical simulation are shown by 
black points, whereas the theoretical values, see Eqs. (4.35) and (4.36), are shown by 



solid line. The numerical values are obtained by averaging over 200 realizations, each 
of length 1024. The vertical marks indicate the velocity relaxation time r^ = Xjav and 
coordinate relaxation time t^ = vJoluj^ for the overdamped oscillator, correspondingly, 
as well as relaxation time t^ = t^ = Ijotv for the weakly damped oscillator. Both 
theoretical and numerical curves reach "plateau" at the intervals greater than the re- 
laxation times. It implies that the processes vif) and x(fy become stationary ones. The 
figures demonstrate an important difference between overdamped and weakly damped 
oscillators: for the overdamped oscillator coordinate relaxation time is much greater 
than the velocity relaxation time, whereas for the weakly damped oscillator both times 
are the same. From Fig. 11 we also make a conclusion about good quantitative agree- 
ment between theoretical results obtained by solving FSFPE and numerical solution of 
the corresponding Langevin equations. In this respect it is worthwhile to point ones at- 
tention to coincidence of theoretical curve and numerical dependence at nonstationary 
parts. 

In order to estimate the influence of power asymptotics on the evolution of the 
moments, we replace stable PDFs in the Langevin equation by "truncated" ones, that 
is those, in which the large values of random quantities are cutted off. In Fig. 12 
the velocity moments versus time are depicted for the weakly damped oscillator in a 
log - log scale. The oscillator frequency is 0.05, the friction coefficient is 0.01, the 
order of the moment is 0.25, and the Levy index is 1.3. The moments are obtained by 
averaging over 1500 realizations, each of length 2048, thus the total number of points 
is 3 ■ 10®. The mode of maximum value (that is, the most probable value) is of the order 
N^/"-. In the figure the moments obtained numerically are shown by black points. The 
Langevin source A{t) is modelled as a consequence of independent random variables 
possessing truncated stable PDF, that is \A{t)\ < A^ax = 1600; 300; 50 for the cases a; 
b; c, respectively. The solid line indicates the moments estimated analytically from the 
FSFPE. It is seen, that the role of large outliers increases with time increasing, since 
large values become more and more probable. Therefore, the discrepancy between 
theoretical results and numerical simulations using truncated PDF grows with time 
growing. With the truncation parameter decreasing the discrepancy increases. Thus, 
as it is clearly seen, the discrepancies are most essential at the stationary stage of 
evolution. 

We have already mentioned that our studies demonstrate a good quantitative 
agreement between theory based on FSFPE and numerical simulations based on the 
Langevin equations. In order to show this fact more precisely, we estimated velocity 
and coordinate moments by averaging over 50-10^ realizations with the total number 
of points 10^, which is much larger than in simulations presented in Fig.ll. Here we 
use the following parameters: u = 0.05, z/ = 0.01, g = 0.25, a = 1.3. By black points 
the results of simulation are shown, whereas the solid lines indicate theoretical values 
estimated with using Eqs. (4.35) and (4.36). The theoretical values of r^, and 5t^ 
are indicated by arrows. It is seen that numerical results strictly repeat all bends of 
theoretical curves at the non -stationary stage of evolution. 

6 Results 

A large section of statistical physics deals with evolution of the systems infiuenced 
by random Gaussian forces. In this paper we study linear relaxation of the systems 
infiuenced by random forces which are distributed with symmetric stable probability 



laws. Stable laws (as the Gaussian one) are the limit ones for probability laws for 
sums of independent identically distributed random variables. Therefore, they appear 
in problems, whose result is determined by the sum of a great number of independent 
identical factors. 

The main results are as follows. 

1. We get fractional symmetric Fokker - Planck equation (FSFPE) and fractional 
symmetric Einstein - Smoluchowski equation (FSESE). They generalize the Fokker - 
Planck and Einstein - Smoluchowski equations for the Brownian motion. The FSFPE 
describes a linear relaxation in the phase space of the systems influenced by stochastic 
forces distributed with symmetric stable law. The FSFPE contains fractional velocity 
derivative instead of the second one. The FSESE describes relaxation in a real space. 
It contains fractional space derivative. To get these equations we use an approach 
analogous to that used by Chandrasekhar in order to derive the kinetic equations for 
the Brownian motion. 

We restrict ourselves by the case of one - dimensional symmetric stable PDFs. The 
approach can be generalized in order to get the kinetic equations for two - and three - 
dimensional real space as well as for the case of asymmetric stable PDFs of stochastic 
external forces. At these ways the problems appear, which are connected with, at 
first, many - dimensional generalization of stable probability laws and, at the second, 
with many - dimensional generalization of fractional derivatives. On the other hand, 
the case of asymmetric stable PDFs may appear to be important for the description 
of asymmetric diffusion. These problems in relation with the kinetic equations are 
discussed in some recent papers [^ , [^ . 



2. With the help of kinetic equations obtained we consider the processes of a linear 
relaxation for two problems: force - free relaxation and relaxation of a linear oscillator. 
We get general analytic solutions of FSFPE and FSESE, and expressions for fractional 
velocity and coordinate moments as well. 

3. For the both problems we solve numerically Langevin equations with the random 
source, which is a discrete approximation to a white Levy noise. After averaging over 
many realizations we estimate fractional moments and compare numerical results with 
the results of analytical solutions to the kinetic equations. Analytical and numerical 
results appear to be in a quantitative agreement. 

4. The process of a linear relaxation in a force - free case can be divided into 
two stages which are described in the framework of FSFPE: "fast" stage, at which 
stationary stable PDF over velocity is established, and "slow" diffusion stage, at which 
relaxation in a real space occurs. The latter process can be described asymptotically 
as the Levy stable process with independent increments. The characteristic time of 
velocity relaxation is r^, = 1/az/, where v is the friction coefficient in the Langevin 
equation, and a is the Levy index of stable PDF of the random force in the Langevin 
equation. The characteristic time of relaxation in a real space is t^ = (z/L)"/D, where 
L is an external size of the system considered, and D^^" is the scale parameter of the 
PDF of the stochastic force. For the large enough systems r^, >> r^, and division 
of relaxation process into two stages is valid. For this case at the diffusion stage 
relaxation can be described in the framework of FSESE, which describes anomalous 
super diffusion, that is, the diffusion process when the particle displacement grows as 
t^/", < a < 2. 

5. When studying relaxation of a linear oscillator, it is expedient to distinguish 
between two cases: overdamped oscillator, w/z/ << 1, and weakly damped oscillator, 
uj/u » 1. Both cases are of special importance in nonlinear generalizations of the 



theory presented. We study in detail, both analytically and numerically, both limit 
cases, and point attention to substantially different properties of relaxation processes 
in two cases. 

6. The relaxation of an overdamped oscillator occurs in two stages, which are 
described in the framework of FSFPE: "fast" stage, at which during time interval 
r„ = l/au stationary stable PDF over velocity is established, and "slow" diffusion 
stage, at which during the time interval r^ = vjauJ^ stationary stable PDF in a real 
space is established. At the diffusion stage the relaxation of an overdamped oscillator 
can be described in the framework of FSESE. 

7. It is known that for a weakly damped Brownian oscillator a method occurs, 
which allows one to simplify kinetic description by passing to slowly varying (at the 
period of oscillations) random variables. We generalize this approach to the case of 
a weakly damped Levy oscillator and get kinetic equation for the PDF depending on 
slow variables. This equation generalizes the kinetic equation for a weakly damped 
Brownian oscillator and contains fractional derivatives over velocity and coordinate. 
Its structure is more simple than that of FSFPE, therefore, it can be applied for non 
- linear problems. We find stationary solution of the obtained equation and show that 
relaxation process can not be divided into two steps. Both velocity and coordinate 
relax during the same time interval r = Ijav. 
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FIGURE CAPTIONS. 

Fig.l. C^/'^{q] a) versus q for different a = 1.0; 1.3; 1.6; 1.9. 

Fig. 2. Force - free relaxation in the framework of FSESE. The moment M^ versus t 
at different values of moment exponent q. The Levy index a is 1. 

Fig. 3. Force - free relaxation in the framework of FSESE. The exponent /i of the 
time dependence of the moment M^, see Eq. (3.13), versus the Levy index a. The 
moment order is g = 0.25. Theoretical dependence 0.25/a is depicted by dotted line. 
Numerical results are depicted by black points. 

Fig. 4. Relaxation of a linear oscillator in the framework of FSESE. The g - th 
order coordinate moment versus time in a log - log scale. The results of numerical 
simulations are depicted by black points, the moment obtained from FSESE is shown 
by solid line. Vertical marks indicate the time of coordinate relaxation. 

Fig. 5. Force - free relaxation in the framework of FSFPE. The typical trajectories 
of the velocity (at the right) and of the coordinate (at the left) for various values of 
the Levy index a. The friction coefficient z/ is 0.03. 

Fig. 6. Force - free relaxation in the framework of FSFPE. The typical trajectories 
of the velocity (at the left) and of the coordinate (at the right) for various values of 
the friction coefficient u. The Levy index is 1.3. 

Fig. 7. Force - free relaxation in the framework of FSFPE. The velocity moments 
(at the left) and the coordinate moments (at the right) versus time in a log - log scale. 
Black points indicate results of numerical simulation of the Langevin equations, the 
solid lines indicates the moments obtained from FSFPE. The moment order q is 0.25, 
The Levy index a is 1.3. Horizontal dots show the time r„ of the velocity relaxation 
and the time 5r„. The dotted lines at the right figures indicate theoretical values q/a 
of straight - line asymptotics (in a log - log scale) for the moment of the coordinate. 

Fig. 8. Relaxation of linear oscillator in the framework of FSFPE. Typical trajecto- 
ries of the velocity (top) and the coordinate (bottom) for overdamped (left) and weakly 
damped (right) oscillators. The frequencies are 0.003 (overdamped) and 0.3 (weakly 
damped) correspondingly, the friction coefficient is 0.03, the Levy index is 1.3. 

Fig. 9. Trajectories on the phase plane for (a) weakly damped, and (b) overdamped 
oscillators, respectively. The parameters are u = 0.06 and uj = 0.009 for (a) and (b), 
respectively, z/ = 0.02, and a = 1.2. 

Fig. 10. Trajectories on the plane [x, E), E is the oscillator energy, for two Levy 
indexes a = 1.9 (above) and a = 1.6 (below), respectively. The parameters are u = 
0.003, u = 0.002. 

Fig. 11. Relaxation of linear oscillator in the framework of FSFPE. The velocity 
moments (at the top) and the coordinate moments (at the bottom) versus time in a log 
- log scale. At the left the results for the overdamped oscillator are presented, uj = 0.01, 
whereas at the right the results for the weakly damped oscillator are presented, uj = 0.1. 
The frequency coefficient u is 0.03, the order q of the moments is 0.25, the Levy index a 
is 1.3. Black points indicate results of numerical simulation of the Langevin equations, 
the solid lines indicates the moments obtained from FSFPE, see Eqs. (4.35) and (4.36). 
The arrows show velocity and coordinate relaxation times. 

Fig. 12. Relaxation of linear oscillator in the framework of FSFPE. The velocity 
moments versus time in a log - log scale for different truncation parameters Amax for 
the stable PDFs of the noise in the Langevin equations. The parameters used in 
simulations are as follows: uj = 0.05, z/ = 0.01, g = 0.25, a = 1.3. Amax is 1.500, 300 
and 50 for cases a,b,c, respectively. 



Fig. 13. Relaxation of linear oscillator in the framework of FSFPE. The velocity 
moments (top) and coordinate moments (bottom) versus t in a log - log scale. The 
parameters of the simulations are as follows: oo = 0.05, z/ = 0.01, q = 0.25, a = 1.3. The 
moments are estimated by averaging over 50-10^ realizations each of the length 2048. 
The vertical marks indicate relaxation times r^,, r^ and 5r„, 5r^, as well. 
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